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Abstract 

It has recently been discovered that for certain rates of mode-exchange cohisions analytic 
solutions can be found for a Hamiltonian describing the two-mode Bose-Einstein condensate. 
We proceed to study the behaviour of the system using perturbation theory if the coupling 
constants only approximately match these parameter constraints. We find that the model is 
robust to such perturbations. We study the effects of degeneracy on the perturbations and find 
that the induced changes differ greatly from the non-degenerate case. We also model inelastic 
collisions that result in particle loss or condensate decay as external perturbations and use this 
formalism to examine the effects of three-body recombination and background collisions. 



1 Introduction 

A Bose-Einstein condensate (BEC) is a state of matter in which a large number of bosons occupy the 
same quantum mechanical ground state. As such, BECs present the opportunity to study quantum 
systems which display large scale (macroscopic) collective behaviour. Recently there has been 
much interest in multi-component BECs because of their importance to quantum optics [H El E], H] . 

*Publishcd before under maiden name Fuentes-Guridi. 



Multi-component BECs are most often formed in a multi-well potential in which the components 
are spatially separated [21 H] . Alternatively, the multi-component formalism can be used to model a 
single component BEC that possesses several internal degrees of freedom, such as varying amounts 
of spin [3] . 

In general, many-body systems are of significant importance in physics. For example, quantum 
information processes require the manipulation and control of systems with large numbers of par- 
ticles. However, many-body systems are often difficult to treat exactly and are most often studied 
using numerical or approximate methods. The realm of applicability of these methods is limited by 
the number of degrees of freedom of the system. Two common approximate models are the Bose- 
Hubbard model of quantum optics [1] and closely related to it, the Lipkin-Meshkov-Glick model of 
nuclear physics [5], both used to describe the two-body interactions of spin- J systems. 

A family of exactly solvable many-body systems was introduced in |6] and studied in greater 
depth in [7j. These models can be used to describe the physics of a two component BEC where both 
elastic and mode-exchange inelastic collisions are present. Mode-exchange collisions are known as 
general nearest neighbor interactions in the context of multi-well BECs [9] and as inelastic collisions 
in the case of a single condensate consisting of particles in two hyperfine levels ^U\. The family of 
models is parameterized by a positive integer n, and hence the specific models are called n-models. 
This terminology is used because the n-model contains 1, 2, . . . , n— body interactions. Analytic 
solutions can be found for the n-model when the strengths of the various particle interactions obey 
specific constraints. While the n-models are of interest to many-body physics in general, in this 
paper we will largely be concerned with the 2-model where only single body interactions and two- 
body collisions are considered. Microscopic calculations show that mode-exchange collisions occur 
in BECs as a result of the interaction of a laser field with the system [10]. The 2-model includes 
the usual Josephson-type interactions, but also allows the effects of mode-exchange collisions to be 
studied analytically, and as such provides a more realistic framework for studying two-mode BECs 
than the canonical Josephson Hamiltonian [n| . 

The analytic solution found in [6j requires that the strengths of the various interactions meet cer- 
tain constraints; if these constraints are not met the solution is invalid. While in some experimental 
situations the rates of inelastic and elastic mode-exchange collisions can be controlled externally 
[12] , this is not always the case. Moreover, even in if the rate of collisions can be manipulated, 
the constraints will likely still only be approximately satisfied. It is thus of interest to extend the 
solution space of the system to the case in which these constraints are only approximately satisfied. 
This naturally leads to the consideration of small parameter perturbations in the model to study its 
robustness. Along with perturbations within the model it is of interest to include additional interac- 
tion terms as perturbations, such as inelastic collisions resulting in particle loss. Particle loss, often 
suppressed in experimental settings, is studied theoretically by considering classical rate equations 
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[TT] . More recently, a quantum treatment of 3-body loss has been given in [13]. Including external 
perturbations extends the predictive power of the model. Additional interactions of primary inter- 
est include three-body recombination, background collisions, spin exchange and dipole relaxation 
[TT] . To our knowledge, this is the first study which analyzes particle loss as a perturbative effect 
in BECs. 

Motivated by the above considerations, in this paper we carry out a full perturbative analysis 
of the solvable BEC 2- model proposed in [6|. We begin by reviewing this model and the solutions 
derived when the constraint equations on the parameters are satisfied [7|. We then introduce 
perturbations to the parameters in the model and analyze the effects on the particle distribution, 
entanglement and the evolution of the relative number operator. In later sections we discuss the 
effects of state degeneracy and the inclusion of a general loss term as a perturbation, illustrating 
the formalism by studying background collisions and three-body recombination. 

2 A Model for Two-Mode Bose-Einstein Condensates with 
Mode-exchange CoUisions 

The 2-model studied in [6, 7] is governed by the Hamiltonian 

H2 = Ao + iu {a^a - b'^b) + X {e"''a^b + e-'^ab'^) + Ua^ab^b 

(1) 

+A [e^"t'a)a)bb + h.c) + /i {e*"^ [a)a)ab - a^b^bb) + h.c] 

where h.c. denotes the Hermitian conjugate of the preceding term. The two modes a and b are 
independent Bose operators satisfying [a, a^] = 1 = [b, b^] with the commutators of all remaining 
pairs vanishing. The term u (a^a — b^b) is the free energy of fia = a^a particles in the a mode 
and fif, = b'^b particles in the b mode, with frequency difference uj between modes. This frequency 
difference arises because the model describes atoms in different hyperfine levels, or alternatively 
unsymmetric spatially separated condensates p3J. A Josephson-type or spin flip interaction is 
included with strength A and phase 0. In practice such an interaction is induced by an external 
field, such as a laser [8]. This interaction may also be interpreted as modeling the tunneling of 
particles between modes with probability proportional to A. Also included are terms corresponding 
to number-preserving elastic and mode-exchange collisions, namely those terms containing four 
bose operators. The interaction with strength yU is a single dispersive process. The interaction with 
strength A describes a collision where two particles exchange their mode. The elastic interaction has 
strength U and models the collision of a particle from each mode in which the number of particles 
in each mode is conserved. We see that the mode-exchange collisions preserve the total particle 
number but not the relative particle number. 

The Hamiltonian given in eq. ([T]) can be efficaciously studied by first introducing the simpler 
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Hamiltonian [6] 



Ho := Ho^2 = Ai {a^a - b^b) + A2 {a^a - b^b)' 



(2) 



with real constants Ai and A2. Such a Hamiltonian models a two-mode condensate with energy 
difference Ai between modes and elastic scattering probability proportional to A2. As is clear from 
the form of Hq the relative number operator m = a^a — b^b is a commuting observable so that the 
number of particles in each mode is conserved. In particular, there is no probability of spin-flip or 
tunneling between modes. We also note that the total number operator = a^a + b^b commutes 
with Hq. We may thus take the eigenvalues and m of N and m as labels of the eigenstates, for 
which we write \N,m). For fixed N E N the values of m are limited tcl] — A^, — A^ + 2, . . . , N — 2 
and A^. The energy of the state | A^, m) is Em = Aim + A2m?'. 

For a certain choice of parameters in H2 the solutions | A^, m) of can be used to obtain analytic 
solutions of H2. To see this, define a two-mode displacement operator by f/ (^) = exp {^a% — ^*ab^) 
with displacement parameter = ^9e^'^. It is clear that U is unitary. It can then be shown that 
if the parameters in H2 satisfy 



Ao = A2 {N^ cos^e + ATsin^ 6) 

u = Ai cos 9 
X = Ai sin 6 
U = 2A2 (1 -Scos^^) 
A = A2 sin^ e 
yU = 2A2 COS 9 sin 9 
then H2 = U'^HqU, a result shown by computing 

/ 
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Observe that U a b at 6t y f/t 

may be computed by setting 9 1— —9 in eq. (jlj). 

Since the eigenvectors of Hq are of the form \N,m) the eigenvectors of H2 when satisfying eqs. 
(|3]) are simply W\N,m) with energy Em = Aim + Aim^. An extensive analysis and discussion of 
these solutions is included in [7]. 

For fixed A^ G N the ground state of 7^2, labeled as \ipmo) = t^t|A^, mo), is found by minimizing 
the energy Em = Aim + A2m^ with respect to m. First assume A2 < 0. In this case, if < 
[Ai > 0) the minimum occurs at m^ = N (mo = — A^). Secondly, say A2 > 0. If 



2A2 



^See the appendix for details. 
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then mo is the closest allowabl^^ integer to — otherwise mo is the closest allowable integer to 
— sgn(y4i)A^. In particular, a choice of either mo or ^ determines the otherjf 

It is important to recall that, in the case of a double-well BEC, the two-mode approximation 
must be satisfied p^. This means that collisions taking place in the region where the wavefunctions 
overlap must be less probable than collisions between particles belonging to the same well. It is 
possible to find an exact analytic solution to the model in this case by considering 6 << 1. It is 
interesting to observe that it is not possible to find analytical solutions to this model for an exactly 
symmetric double-well and satisfy, simultaneously, the two mode approximation. Fortunately, in 
the case of slightly asymmetric double wells (which corresponds to a more realistic situation) it is 
possible to find exact analytical solutions in the two-mode approximation [7] . 

The most general n-model Hamiltonian discussed in |6] has as its Hamiltonian Hn = U^Ho^nU 
where ifo,n = Sfc=i^fc^'^- The eigenvectors of Hn are again f/^|A^, m) with the same energy Em- 
The same process of matching coefficients in the expansion of WHo^nU that was used to obtain 
eqs. ([2]) can be carried out so that W\N,m) is the exact solution of a Hamiltonian containing up 
to n-hodj interactions. In the general case, there are n + 2 free variables (the n Ai as well as 6 and 
0) in UHo^nU^ , whereas the number of terms in a general Hamiltonian describing such interactions 
is (n + 2)(n + 4)(n + 6)/48 as shown in Appendix B. So while the above method continues to give 
analytic solutions for all n, its range of applicability decreases as n grows. It should also be noted 
that because Ho n contains only terms with the same number of creation and annihilation Bose 
operators, this formalism cannot be used to study interactions involving an odd number of Bose 
operators, such as the interaction a^ba where one particle is lost. 

We remark that all of the perturbative analysis carried out below is valid not only for the 2- 
model, but also for the more general n-models. This is so because the solutions for all n-models are 
U^N,m). We need only change the interpretation of the perturbative analysis in each case. 

3 Effects of Parameter Perturbations in the 2-Model 

While satisfying eqs. is sufficient to obtain an analytic solution of H2, the method presented 
above fails to produce solutions if the parameters in H2 deviate even slightly from these constraints. 
Thus, in order to compare this model to experimental results we would like to study the solutions 
to H2 when the constraints are only approximately satisfied. We proceed with these perturbation 
calculations in the section below. 

We begin by perturbing each of the coupling constants in H2 away from the conditions given by 
eqs. ([3]), i.e. perturbations of the form u u + 6uj where u satisfies eqs. ([3]) and 6^^ is small. We 
omit the study of perturbations Ao Aq + 6ao since this simply leads to a shift in the energy of the 



^The integer m can only take the values —N, ~N + 2, . . . , N — 2 and N; see the appendix for details. 
•^In calculations we will assume A2 = 1, so that a choice of toq determines Ai. 



5 



state by 6ao- In what follows, we assume that all eigenstates are non- degenerate; the degenerate 
case is discussed in Section 7. Throughout the paper, given an operator O we define O := UOW. 
To perform many of the calculations that follow we have made use of a coordinate transformation 
between the {N,m} basis and the {ua^Ub} basis, \N,m) i-^ \na — ^^^^,^6 = ^^^^^) where is a 
eigenvalue of fia — a'a and similarly for rib- 

3.1 LJ Perturbation 

A change a; i— > a; + 5^; results in the perturbation H^^ — 5^ [a^a — b^b) . A calculation then shows 
that the non-vanishing matrix elements of are 

{N,m\H^\N,m) — 6^m cos 6 (5) 

and 

r 

{N, m ± 2\H^\N, m) = -^e^'^ sin OyJ N{N + 2) - m{m ±2). (6) 

3.2 A Perturbation 

A change A i— > A + 5^ results in the perturbation E.\ — bx {e^'^a^b + e"~"^a6^) and we find that 

{N,m\Hx\N,m) = 5xm sin 6 (7) 



and 



{N, m±2\Hx\N,m) = ^ ( e^*^ cos^ ^-6 - e^^'^ sin^ J ^J N{N + 2) - m{m ±2). (8) 



3.3 U Perturbation 

A change U + 5u results in the perturbation Hu — Sua^ab^b and we find that 

{N, m\Hu\N, m) = I sin^ 9 (^^^l±Ill + + Su cos' 9 (^^^^-^ ] , (9) 

{N,m±2\Hu\N,m) = -^e^*^ sin^ 9y/{N ±m) {N ±m + 2) {m ±1) (10) 

8 

and 

(N,m±i\Hu\N,m) ^ — -e^^''^ sin^ 9y/[N(N + 2) - m(m ± 2)1 [N(N + 2) - (m ± 2)(m ± 4)1. 

16 

(11) 
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3.4 A Perturbation 

A change A i— > A + 5a results in the perturbation = 6^ (e^*'^a"l"a'^66 + h.c.) and we find that 

{m\HA\m) = I sin^ 6 ^t!^^ _ , (12) 

{m±2\Hj,\m) = ^e^'^ sin 29 ^/{N t m) {N ± m + 2) {m ±1) , (13) 

and 

(m ± 4|ifA|m) = ^e^^*"^ (l + cos^ 6) ^/{N ±m + 2) {N ±m + A) {N t m) {N t m - 2). (14) 
8 

3.5 fi Perturbation 

A change jj, ^ jjL + 5^ results in the perturbation = 5^ {e*'^ {a^a^ah — a'^b'^bb) + /i.e.} and we find 
that 

{m\H^\m) = ^sm29i Nj, (15) 

r 

±2\H^\m) = -^e^'"^ cos 29 ^/{NTm) {N ± m + 2) {m ± 1) (16) 



{m 
and 



{m±A\H^\m) = — i^e^^'^ sin 29^y{N ±m + 2) {N ± m + A) (N t m) (N t m - 2). (17) 
8 

In each of the cases above the perturbation matrix elements simplify significantly in the limiting 
cases 9 = 0,71, in which the Josephson coupling vanishes (A = 0), and 6* = | in which the potential 
well is symmetric {uj = 0). We elaborate below on the effect of 9 on the particle distributions. 

4 Perturbative Effects on Particle Distribution 

In the unperturbed case we can use the analytic solution \ipmo) = U^N, uiq) for parameters satisfying 
eqs. (13]) to find an explicit expression for the particle distribution: 

P =\{N,m\U^\N,mo)\'^. (18) 

Using the homomorphism described in the appendix to relate the Schwinger su(2) representation 
to the angular momentum representation we write P = {d^j^J"^ wherj^ [15] 



<m. = (14.], l^j), [-Os'-] S (19) 



mo '"^^^ Wigncr rotation matrix elements under the effect of the Lie algebra homomorphism discussed in 

the appendix. To see why these appear, observe that C/t = g^*^" '' where n — (sin0, cos 0) and J = (^Jx, Jy, Jz^ 

is the total angular momentum vector operator. That is, U\ and hence U, are rotations of the algebra generated by 
a and b. 
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Figure 1: The figures show the unperturbed particle distributions P given by eq. (fTSi) for = 1000, 
^ = 1 and a) mo = 1000, b) mo = 998 and mo = 996. 
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(20) 



The integers k± are chosen so that the arguments of the combinatorial symbols are non-negative; 
exphcitly, k. = max {O, ^no^] and k+ = min {^, ^^}. We plot in Figured the unperturbed 
particle distribution for = 1000, 6 = 1 and mo = 1000, 998 and 996. Observe that the particle 
distribution is independent of phase (p. We will see that this does not hold in the perturbed case. 
As 6 varies the distributions shift along the m-axes; for 6 small the maxima shift toward m = N, 
and move toward m = — A^ as 6 grows. 

The canonical 2-mode BEG predicts (under certain circumstances) that the ground state solution 
is a superposition of two peaked distributions. From Figure [1] we see that for mo < N — 2 the ground 
state is a superposition of more than two distributions, an effect that is due to the mode-exchange 
collisions not included in the canonical model [7]. 

In order to include perturbative effects in the particle distribution we replace the zeroth order 
ground state wave function W\N,mQ) in eq. f[T5]) with that including first order corrections, 

f/t|A^,mo)(°+^) = {\N,mo) + |A^,mo)(^)) . 



Recall that in the non-degenerate case a perturbation H' induces a first order wave function cor- 
rection given by 



\N,mo) 



(1) 



V- . V- (Ar,m|g^|Ar,mo) |,^ 
2^ a^o,^|A^,m) := 2^ |A^,m) 



E - E 



(21) 



where the superscripts indicate to what order in the perturbation the term corresponds; when no 
superscript is present the result is to zeroth order. With this, we find that to first order the particle 
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distribution is given bja 

N 

P*- ^ = Ifi^jmol "I" 2*^771, mo ^ ^ -R-^ ('^mo,n) C^m.n- (22) 

n=-Af 

We show in Figures [2] - [6] the particle distributions under each of the five parameter perturbations 
considered above. Note that in certain cases p(°+^) may be negative because the first order term 
in eq. (l22ll depends hnearly on the perturbation coefficients am,n] this of course is just a reffection 
that only first order corrections have been taken into account. 

We see from eq. (1221) and the fact that am,n is proportional to the sign of the perturbation that 
each of the perturbations can be used to either enhance or diminish the particle distributions in 
the vertical direction for a given m. That is, given a fixed relative number m, by choosing the 
sign of the perturbation appropriately we can increase or decrease the probability of the system 
being in the state \N,m). We also see that as |mo| decreases the number of local maxima of the 
particle distribution increases, and the asymmetry of the perturbations becomes more significant, 
vertically diminishing probability amplitudes on one side of the central maxima while vertically 
stretching the amplitudes on the opposite side. The figures also show that the system is much less 
sensitive to perturbations in u and A than it is to perturbations in the collision coupling constants; 
in the figures the perturbations of uj and A are between 7 and 100 times greater than those in the 
collision perturbations {A, fi and U), while the corrections in all of the cases are of the same scale. 
Moreover, out of the collision parameters the particle distribution is most sensitive to perturbations 
in coupling constants of the mode- exchange collisions (/i and A). It is thus beneficial to maximize 
the matching of the mode-exchange collision constraints ([3]) through redefinition of parameters, 
allowing the perturbation to have a larger effect on the uj, A and U terms. 

The effect on the particle distributions of varying 6 is to shift the position of the central maxima 
on the m-axis: for 6' = |, corresponding to the case in which the Josephson-type interaction is 
maximal (A = Ai) and the energy of each mode is equal [uj = 0), the distribution is centered 
around m = 0, while for 6' = (6' = vr), in which the Josephson-type interaction vanishes (A = 0) 
and the energy difference of the modes is maximal (|ti;| = Ai) it is centered m = N [m = —N). 
In each of the cases, despite being centered around different values of m the perturbations display 
the same qualitative behaviour for any choice of 6. Throughout the paper we have chosen 6 = 1 
in the figures as this lies between the extremes of 6 = and 6 = ^. As 6 varies the effects of 
the perturbations vary in size, but show the same characteristics. For example, as 6 decreases 

N 

^Throughout the rest of the paper, where appropriate, we write with the understanding that we actually 

mean . Summations of this form arise because of the restriction on the eigenvalues of the 

nG{-N-N+2,...,N-2,N} 

relative number operator. 
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Figure 2: The figures show the unperturbed (sohd) and perturbed (dot) particle distributions P 
given by eqs. ([HD and ([22]) for = 1000, 9 = 1, = 15 and a) mo = 1000, b) nio = 998 and c) 
mo = 996. 




1000 




b)mo = 998 
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Figure 3: The figures show the unperturbed (sohd) and perturbed (dot) particle distributions P 
given by eqs. ([18]) and ([22]) for = 1000, 9 = 1, 5x = 15 and a) mo = 1000, b) mo = 998 and c) 
mo = 996. 

from 1 to (but with 6 fixed), the perturbations have a larger effect on the particle distributions, 
which indicates that for small values of 9, we must be more careful about the size of the chosen 6. 
This can be explained as follows. In the perturbed case the area under the particle distributions 
sums to unity only to first order in the perturbation. This sum is also a function of 9 (since the 
perturbed terms are), unlike in the unperturbed case. We have found that as 9 decreases the sum 
of probabilities decreases too, which explains why there is a difference in the particle distributions 
when we keep the size of the perturbations 6 fixed. If we decrease both 9 and 6, then we find that 
the particle distributions essentially only shift along the m-axis, but not in their general shape. 



5 Perturbative Effects on Entanglement 

Entanglement arises in many-body quantum systems because of the superposition principle and the 
tensor product structure of the Hilbert spact^, a property of utmost importance in quantum control 
and quantum information. In particular, states with high entanglement are desirable because of 
their utility in carrying out quantum information tasks [TB]. We would like to determine whether 

N-m\ 



^Note that we have been using the shorthand \N,m) = \na = ^^"' ) 



\nb 
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Figure 4: The figures show the unperturbed (sohd) and perturbed (dot) particle distributions P 
given by eqs. ([HD and ([22]) for = 1000, 9 = 1, Su = 2 and a) mo = 1000, b) mo = 998 and c) 
mo = 996. 




r/f 



Figure 5: The figures show the unperturbed (sohd) and perturbed (dot) particle distributions P 
given by eqs. and ([22]) for = 1000, 9 = 1, 5a = 0.1 and a) mo = 1000, b) mo = 998 and c) 
mo = 996. 



or not by careful choice of the perturbations we can increase the entanglement from that of the 
unperturbed case. 

For a bipartite quantum system the von- Neumann entropy S = —Tr (p logg p) is a standard 
measure of the entanglement of the system, p being the reduced density matrix. In the case at hand 
this reduces to 

N 

S{N,mo,9) = - P{N,mo,m,9)log,P{N,mo,m,9) (23) 

ni=-N 

where P is either the perturbed or unperturbed particle distribution, depending on the situation. 
As noted in the section above, the first order particle distribution given by eq. fl22]) can be negative 
for certain choices of perturbations. To resolve this problem we could either replace p(°+^) with 

2 

|p(o+i) I Qj. could include the second order term J2n=-N '^rno,nd^^n p(o+i); we use the former 
approach to keep the calculations strictly to first order. 

When the particle distribution contains first order corrections, p(°+^) = p(°) + P^^\ the entan- 
glement is given by 



^ / 1 \ 

^(o+i)^5(o)_ J2P'^^\N,mo,m,9)nog^P^''\N,mo,m,9) + — ]. (24) 

m=-N ^ ^ 
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Figure 6: The figures show the unperturbed (sohd) and perturbed (dot) particle distributions P 
given by eqs. ([18]) and ([22]) for = 1000, 9 = 1, 5^ = 0.1 and a) mo = 1000, b) tuq = 998 and c) 
uio = 996. 



f . A 




Figure 7: The figure shows the unperturbed entanglement using eq. ( l23l) for = 100. 

For the case at hand we read off that P^^^ = 2(i^ Re (amo,n) c^m,n- Hence we have a 

criterion for increasing the entanglement of the system. Namely, we increase the entanglement (to 
first order) precisely when 

P^^KN,mo,m,e) [\og^P^''\N,m,,m,e) + — \ < 0. (25) 

m=-N ^ ^ 

Note that P^^\N, mo, m, 6) is proportional to the perturbation strengths, so that we can increase or 
decrease the entanglement by choosing the sign of the perturbation appropriately. Continuing, we 
note that for each m we have logg P^'^^'^N, mo, m,9) + ■^^ > 0, so that to maximize the entanglement 
we should maximize each P^^\ Of course, |-P^^^| can be made arbitrarily large simply by choosing 
S arbitrary large. However, we are limited in such a choice since we must keep 6 small so that 
perturbation theory can be trusted. 

We plot in Figure [7] the unperturbed entanglement as a function of 6 and mo. See [7] for an 
extensive analysis. In Figures ISlfT^ we plot the perturbed entanglements 5''^°+^^ for 6 = 0.01 as well 
as the differences AS" := S^^~^^^ — S^^^ between the perturbed and unperturbed entanglements for S = 
0.1; we choose a different S in the latter case so that the differences are more evident. In each of the 
plots of AS* we see that the largest changes occur approximately along the diagonal lines, connecting 
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Figure 8: The figures show the a) perturbed entanglement for 6^ = 0.01 as a function of mo and 9 
for iV = 50 and 6; the difference - for N = 100 and 5^ = 0.1. 

the points {6, itlq) = (vr, — A^) and (vr, A^) and the points (— vr, — A^) and (vr, A^); these hues in 6' — mo 
space correspond to certain strengths of the couphng constants in the Hamiltonian H2 viewed as 
functions of rriQ. Note this type of behaviour also occurs along these lines in the unperturbed 
entanglement plots, as shown in Figure [71 Comparing the magnitude of the perturbative effects 
we observe that perturbations to the mode-exchange collision terms {U, A, fi) have a much greater 
effect, their maximum difference being about an order of magnitude larger than those for A and 
uj. It is also interesting to observe that the coherent states, which correspond to U^N,N) and 
U^N, —N) are the states which are maximally affected by perturbations. In such states Ai is much 
larger than A2 so that the rate of collisions is relatively small. To further study this we have plotted 
in Figure [12] a two-dimensional cut of Figures ^ )-[T^ ), where we have fixed = j and allowed mo 
to vary. We observe that in each of the plots AS* has both a local maximum and local minimum as 
mo approaches A^. 

We also observe that the most significant changes occur in the region mo > 0, where more 
particles lie in the a mode. Positive mo implies that, for fixed Ai > 0, the scattering length for 
same-mode collisions is positive and as mo approaches A^, the collision rate becomes smaller. We 
then conclude that condensates with negative scattering lengths are more resilient to parameter 
perturbations and high same-mode collision rates help stabilize the condensate. 

6 Evolution of Relative Population 

With an analytic solution to the system we may study the evolution of the relative population 
(m)(t) := {tp{t)\m\ilj{t)) as a function of time, where the initial state is given by \ilj{t = 0)) = 
Cmt/^|A^, m). For simplicity we restrict our study to the case in which Cm G M. To first 
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Figure 9: The figures show the a) perturbed entanglement for 6\ = 0.01 as a function of mo and 9 
for N = 100 and b) the difference S^°+^^ - S^°^ for N = 100 and Sx = 0.1. 




Figure 10: The figures show the a) perturbed entanglement for 6u = 0.01 as a function of mo and 
e for AT = 100 and b) the difference S^°+'^^ - S^°^ for N = 100 and Su = 0.1. 

order in the perturbation parameter we have 

N f ^ 1 

m=-N t k=-N ) 

We compute (m)(°+^)(i) = (m)(i) + (m)(^)(i) where 

N N-2 

{m){t)^ cose E mC^ + sin^ ^ C'mCn.+s V^(A^ + 2) - m(m + 2)L^(i) (26) 

m=-N m=-N 

with := cos (0 + (-E'm+2 - -E'm) ^) and 

N r N 

m=-N U=-Ar 

(Af-2 
E a^,^C^G+2 + 2) - l{l + 2) cos [0 + (£;,+2 - + (27) 

l=-N 



N 

E am,iCmCi^2^/N{N + 2) - Z(Z - 2) cos [0 - (£;,_2 - E^)t] 

l=-N+2 
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Figure 11: The figures show the a) perturbed entanglement for = 0.01 as a function of mo and 
^ for iV = 100 and b) the difference S^°+^^ - S^°^ for = 100 and 5a = 0.1. 




Figure 12: The figures show the a) perturbed entanglement for 6^ = 0.01 as a function of mo and 
eioT N = 100 and b) the difference - for N = 100 and 5^ = 0.1. 

We plot in Figure [H] a) the unperturbed evolution of relative population given by eq. fl26|) . 
We see the Rabi-like oscillations with relative population collapse and revival. In Figures [H] b ) 
and we plot the evolution of relative population under the parameter perturbations 6^^ = 0.005 
and 5^ = 0.05, respectively. We observe that as 5^ increases the time-averaged value of (m)(°+^) 
decreases; from eq. (P7|) changing the sign of 6^ would have increased this average value. We see 
in Figure [13] c) that the perturbation has broken down; the maximum value of \{m)^^^^^\ is greater 
than the total particle number. This breakdown reminds us that we must not let the perturbation 
grow too large for our analysis to be reliable. We also see (in Figure [T^ c) that as 6^^ grows the time 
of population collapse significantly decreases. Because of the complexity of the correction term eq. 
(!27|) an analytic study of the effects of the perturbations on population collapse and revival times 
is not possible. 

/^From eq. fl26|) we see that in the unperturbed case an initially pure state remains pure. Eq. 
( !27|) shows that the same holds true even in the perturbed cases. Hence, {m)^^~^^\t) has non-trivial 
time-dependence if and only if the initial state is entangled. 
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Figure 13: Shown in a) - e) are plots of 5''-°+^^ — S^^^ for to, A, W, A and /i, respectively, for 6* = | as 
a function of mo- The perturbation strengths used are all S = 0.1. 
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Figure 14: Figure a) plots the unperturbed evolution of relative probability {rh){t) for N = 50. 
Figures &j and c) plot the perturbed evolution {rh)^^~^^\t) given by eqs. fl26l) and fl27|) for = 50 
with = ^ and = ^, respectively. 



7 Degenerate Perturbations 

Throughout the above analysis we have assumed that the unperturbed states in question are non- 
degenerate. We proceed now to study the degenerate case, which results for specific values of the 
constants Ai and A2. For fixed total particle number two distinct states, labeled by relative 
population numbers mi and m2, have the same energy precisely when mi -|- 11x2 = From the 

analysis of the perturbations completed above we know that the perturbation matrix elements are 
non-vanishing only if mi — 1712 G {2,4}, where we have assumed without loss of generality that 
mi > 7712- Combining these two observations we see that there are at most two pairs of degenerate 
states. 

Let us consider as an example the perturbation u u + 6^ with Ai = — [A^ -|- (A^ — 2)] A2, so 
that the only pair of degenerate states is |A^, A^) and |A^, A^ — 2); each have energy y42A^ (2 — A^). 
The matrix of interest is 



{N, N\HJN, N) {N, N\HJN, N - 2) 
{N,N -2\H^\N,N) {N,N -2\H^\N,N -2) 



(28) 



its eigenvalues and eigenvectors yield the first-order energy and wave function corrections. Using the 
calculations performed above it remains to find the solutions e± of the quadratic det (A^; — e/2x2) = 
0. The eigenvalues (energy corrections) and corresponding eigenvectors (wave function corrections) 
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of are 



^ = {N - 1) cose ± ^/N - {N - 1) cos'^ 9 



(29a) 



±) = a± (VNsm9\N,N -2) + ( cos^ ± v/iV^^TiV^^^)^»s2^) |A^,A^) 



(29b) 



with a± a suitable normalization constant. 

Figure fTSl a) plots e± as a function of 9 for = 1000 and S^o = 0.01. We see that the energy is 
lowered for all values of 9; had 6^^ been negative the opposite would have been true. We see from 
Figures [TBI b) and c), which plot the perturbed particle distributions for 6^^ = ±0.01, that even 
for very small perturbations in u there is a noticeable change in the particle distribution. This 
is in contrast to the non-degenerate perturbation of u, where even for 6^ = 15 there was not a 
large change in the particle distribution. The induced change in the particle distribution is also 
qualitatively different from that in the non-degenerate case. This arises because the correction 
in the degenerate case is sinusoidal with frequency much greater than that of the unperturbed 
particle distribution. Regardless of the sign of S^j we see that the central maximum of the particle 
distributions is shifted to smaller values of m. 

We plot in FiguredHlay' the perturbed entanglement as a function of 9 and mo for 6^ = 0.005 and 
N = 100, while b) of the same figure plots the difference AS for the same configuration. Even for 
a small perturbation strength 5^ = 0.005 the perturbation to the entanglement is still significant, 
again showing that a degenerate system is more sensitive to perturbations than is the non-degenerate 
case. We also observe that the largest perturbations are present for mo close to — and 9 ^ j,^. 
For mo close to A^, the perturbations become negligible, regardless of 9. Therefore, we find that 
high collision rates also help stabilize the condensate against perturbations in the degenerate case. 
However, in this case, condensates with positive scattering lengths are more stable. 

Similarly, we find that for a perturbation A t-^ A + 5a, assuming again that |A^, A^) and |A^, A^ — 2) 
are the degenerate states, the energy corrections are 



which is just the energy correction equation for perturbations in uj with the substitution cos 9 
sin^; we obtain the wavefunction corrections |±) for 6x in the same manner. Figure [T71 plots the 
perturbed energy and perturbed particle distributions for 6x = 1 and 6x = ±0.01. The same 
comments made above about Figure [15] for perturbations 6^^ hold in this case as well. In Figure 
[18] we plot the perturbed entanglement and entanglement difference AS for 6x = 0.005. The figure 
shows similar behaviour to Figure [161 However, we note that in the case of A perturbations AS" is 
strictly positive (for positive 6x)- There are again extrema (this time both maxima) for mo = — A^ 
and 9 ^ f , X' wi^^ AS* vanishing as mo approaches A^. Also as in the case above, the system is 
very sensitive to perturbations, with AS ~ 10 at some points. 




'A 



(30) 
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Figure 15: Figure a) plots e+ (solid) and e_ (dot) from eq. fl29al) for degenerate perturbations with 
= 1000 and 6^ = 1 as a function of 9. Figures b) and c) plot the corresponding unperturbed (dot) 
and perturbed (solid) particle distributions P with N = 1000, 9 = 1 for 6^^ = 0.01 and = —0.01, 
respectively. 




Figure 16: The figures show the a) perturbed entanglement for S^^ = 0.005 as a function of mo and 
^ for = 100 and b) the difference ^(°+^) - with 5^ = 0.005 for = 100. 

The study of degenerate perturbations in the remaining parameters is completed in the same 
way as is done above, so we omit these. Note that high collision rates help stabilize the condensate 
against perturbations in A for both positive and negative scattering lengths. 

8 External Perturbations 

As mentioned in the introduction it is of interest to study perturbations that model additional 
interaction terms not included in the original Hamiltonian eq. ([T]). We will largely be interested in 
perturbations that do not preserve the total number of particles in the system. We begin with a 
general discussion of loss terms and proceed to use this formalism to discuss the effects of background 
collisions and three-body recombination. 
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Figure 17: Figure a) plots e+ (solid) and e_ (dot) from eq. (130|) for degenerate perturbations with 
= 1000 and 6\ = 1 as a function of 9. Figures b) and c) plot the corresponding unperturbed (dot) 
and perturbed (solid) particle distributions P with = 1000, 6* = 1 for 5^ = 0.01 and 6x = —0.01, 
respectively. 




Figure 18: The figures show the a) perturbed entanglement for Sx = 0.005 as a function of mo and 
^ for = 100 and b) the difference S^°+^^ - 5^°) with 6x = 0.005 for = 100. 

8.1 General Loss Terms 

Interactions in BECs that do not preserve the total particle number are often minimized in exper- 
imental settings as there is currently no known analytical model that involves such terms. While 
the Hamiltonian H2 we are studying also does not include such terms, in the case that the effects of 
particle loss terms are expected to be minimal, we may treat loss terms as perturbations to the sys- 
tem. A primary source of particle loss is inelastic collisions p!T| [T7]. In magnetic traps particle-type 
exchange terms dominate loss mechanisms, while in optical traps, these terms may be neglected 
[TS] . It is thus of interest to obtain predictions of the effects of loss terms by treating them as 
perturbations to the exactly solvable 2-model. 
The most general loss term can be written as 

ka=0 fc6 = 
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where fka.ki, (^aj"^?)) are some functions and = A^^ + A*";, is the total particle number of the 
condensate. By probabilistic arguments we expect fka,kt ^ as ka + grows, so that the cases 
of primary interest are those with ka and kh small. However, it should be noted that higher order 
collisions (i.e. not just ka = 1,2, /cf, = 1,2) are of physical significance, particularly when the 
condensate is in its coldest phase and of high particle density [TH [IT]. By suitably choosing fkaM 
we can model specific loss terms. 

In order to study loss terms we must first make some adjustments to the analysis performed 
above. The full Hilbert space Ti of the Hamiltonian eq. ([2]) can be orthogonally decomposed as the 
Fock space H = where, in the {A^, m} basis, 7i„ := { \n,m) \ m = —n, —n + 2, . . . , — 2, n}. 

Since the Hamiltonian and perturbations we have considered so far have all commuted with the 
total number operator we have been able to first choose a total particle number A^ for the system, 
or equivalently the subspace T^at C 7i of the total Hilbert space, and then proceed with calcula- 
tions. In order to study loss terms we must enlarge the state space to be 7i_A = ©„g^'^n where 
^ {0}U^ is of accessible total particle numbers. For example, if Hioss oc then 



The energy of the state | A^, m) is Em = Aim+A2m?. Although Em is only functionally dependent 
on m, it has an implicit dependence on A^ since A^ restricts the values of m. Hence, while \N,m) 
is non-degenerate as an element of TYtv, it may be degenerate as an element of 7i_4 depending on 
m and A. Again, considering the example in which Hioss oc we see that \N,m) G Hn © 'Hn-2 
is non-degenerate if m = ±A^ and degenerate otherwise. In general, let iS^(m) = {|n, m) | n G A} 
and say that iS_4(m) is degenerate if it contains more than one element and say it is non-degenerate 
otherwise. We examine the effects of the degeneracy iS^(m) below. Note that the degeneracy 
studied in this section is caused by the interactions (which determine the Hilbert space), whereas 
the degeneracy studied in the previous section was caused by a specific choice of the coupling 
constants Ai and A2. 

In order to study the particle distribution of the perturbed state |A^, mo)*-"^^-* we must modify 
eq. ffTSl) : if we were to use this formula there would be no perturbative effects on the particle 
distribution since Tij and Tij are orthogonal if z 7^ j. A suitable generalization is given by 



Observe that in the case of no loss terms A = {N} so that Pgen reduces to eq. (ITS]) . 
8.2 Effects of the Degeneracy of <S^(m) 

If iS^(m) is degenerate it is easy to see that the matrix of interest to degenerate perturbation theory 



A={N,N-2}. 




2 



(32) 




is triangular with zeros along the diagonal. Indeed, Hioss contains only 
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annihilation terms, and conjugation by U, denoted here by ~, maps annihilation operators to 
annihilation operators. So Hioss\nj, m) is a sum of states with total particle number n less than rij, 
showing the matrix at hand is triangular, and thus has a trivial spectrum consisting of only zeros. 
Hence we can learn nothing from first order perturbation theory. 

Alternatively, if iS^(m) is non-degenerate we may apply the tools of non-degenerate perturbation 
theory. Although it is straightforward to compute the matrix elements of Hioss given by eq. ( l3Ti) 
in general, the requirement that Sj[{m) be non-degenerate severely limits the usefulness of such a 
calculation. We instead focus on some specific choices of fkaM that model interactions of physical 
interest . 

8.2.1 Background Collisions 

Background collisions most often occur in BECs when particles from the condensate collide with a 
residual background gas in the condensate chamber, or alternatively, with metastable atoms within 
the condensate [19] . Background collisions become more important as the density of the condensate 
increases. As a simple illustration of how we may treat background collisions as perturbations, 
consider the case in which the initial state is |A^, A^), so that only a mode particles can be ejected. 



The general loss ternjll eq. ( l3Ti) . after conjugation by U, is written as Hioss = J2k=i ^'^'^ \Qa^ + 
0(6) where we use 0(6) to denote terms with more h powers than fe"!" powers. Note that any term 
of 0(6) annihilates |A^, A^) so that the perturbative correction can be found by neglecting all such 
terms. We compute the desired matrix elements to find that 



J^a,cos'=Wn;io(A^-j) 



Ml + k(2N - k)A2 



which yields the generalized particle distribution 



«fcCOS^|./n-"n (A^-j) 
^Sen - \d^,n,N\ +^2^ kA, + k{2N - k)A, "^^'^"^rn^N -k l^4j 



k=l 



where we set d^~j^_j^ = if |m| > N — k. From this expression we see that we could have omitted 
terms that eject an odd number of a-mode particles. Indeed, if k is odd, then d^ jyd^J^,^ is 
identically zero as a function of m since if c?^ ^ 7^ then A^, m and n all have the same parityo We 
plot Pgen in Figure Uni considering terms that eject 2,4 and 6 particles for A^ = 1000 and ^ = 1. 
In Figure [inia^ we set aa = -0.1, = -0.001 and ag = -5 x 10"^ In Figures [H]?); - we 
increase each of the by a factor of 2. The figures show that increasing each decreases the 
height of the particle distribution. Also evident from the figures and the scale of the ai is that as i 



^Explicitly, we take fk^.kt = if fcf, 7^ and fk^.o = Oika ^ ^ otherwise. 

®We can, however, use second order perturbation theory. In this case terms annihilating an odd number of 
particles will have an effect on the particle distribution. 
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Figure 19: The unperturbed (solid) and perturbed (dotted) particle distributions, the latter given 
by eq. ([31]), for N = 1000 with a) a2 = -jq, = -j^ and = -2o^> ^y* "2 = -jq, 

'^4 = ~1000' '^6 = ~ 200000 ' '^2 = —10) ^^4 = ~iooo '^6 = — 200000 ^-^ '^2 = —Jq, ^4 = ~iooo 

and tte = — 200000 ■ 




Figure 20: The figures show the a) perturbed entanglement for 0^2 = — ^, q;4 = and = 

— 200000 (as in Figure [T9l a ) as a function of thq and 6^ for = 100 and b) the difference 5*^°"*"^^ — S^^^ 
for the same configuration. 

increases the perturbations have a larger effect. That the perturbations do not blow up reflects the 
requirement that fka,kt ~^ for large ka + kfy. Indeed, we have found using numerical simulations 
that the perturbations that eject k particles diverge with increasing k. 

We remark here that as a consequence of d^ j^d^-^_^. vanishing for k odd we cannot learn 
anything about spin-flip terms from flrst order perturbation theory. Indeed, such terms would be 
modeled by perturbations of the form li^aa and a^aa, assuming particles in the a mode have greatest 
energy; these terms clearly reduce the total number of the system by an odd number. 

Setting P(^) equal to the second term on the right hand side of eq. flM|) we can use eqs. f l24] - 
[251) to compute and increase the entanglement, respectively. Since P^^-* depends linearly on the 
interaction strengths Uk we see the most obvious manner in which to increase the entanglement is 
to make \ak\ large; the sign of ak depends on fc, 6 as well as Ai and A2. 

Figure [20] plots the perturbed entanglement caused by the inclusion of background collisions 
as a function of mo and 9. We see that the largest effect of the perturbation occurs in the region 
6' ~ I and small (large negative) values of mo. There is also a minimum of smaller magnitude 
around 9 = ^ and mo = — A^. The region with large negative mo corresponds the case where the 
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scattering length between particles is negative, in which most particles lie in the b mode of the 
condensate. Since the background collisions considered here eject particles from the a mode, they 
serve to further decrease the value of m. We understand the large effect of the perturbation on the 
aforementioned region as follows: since most particles lie in the b mode, ejecting any particles from 
the a mode has a large effect on the system, since it already has only a small number of a mode 
particles, relative to the number of b mode particles. Figure [20] b) shows that for values of mo above 
the region in question, where the b mode particles become more scarce, the perturbation has little 
effect on the system. Again, we understand this as being because ejecting an a mode particle from 
a state with large m value is of little significance to the system as a whole. We also see from this 
figure that the entanglement is decreased, regardless of mo and 6, for this specific choice of a^. 

8.2.2 Three-Body Recombination 

Perhaps the most physically important type of inelastic collision leading to particle loss is three- 
body recombination (TBR) [20j . TBR occurs when three particles in a single mode collide to form 
a diatomic molecule and a particle of the same mode that carries off any excess energy. Depending 
on this energy and the energy of the potential trap, the resultant particle may or may not escape 
the trap [20] • 

A perturbation modeling TBR would be most naturally treated in the su(3) formalism where the 
extra bose operator would correspond to the diatomic molecule. However, in the su(2) formalism 
we have no such third mode. We thus model TBR by the Hamiltonian 

Htbr = C [aa^aaa + (1 — a)aaaj . (35) 

The parameter a describes the probability of the emitted particle remaining trapped in the conden- 
sate. The term proportional to aaa lowers the total particle number by an odd number and hence, 
as explained above, will have no effect on the perturbed particle distribution. We thus neglect 
this term from our analysis and absorb the constant C into a. With Htbr as a perturbation we 
have A = {N, N — 2, N — 3}. We take as our unperturbed state |A^, A^), which is non-degenerate. 
Proceeding, we find 

^t,. cos^ ^ J N(N -1)(N -2), ^ , ^ 

Ar,Ar«=(T ^-^ — ^ , , -\N-2,N-2) 36 

' ' ' 2Ai + 4:A2{N - 2) I ' / V ^ 



which gives the generalized particle distribution 



Pgeni'm) — \d^^,y\ +2(7 + 4^2 (A^ — 2) ^m,Nd'm,N-2- (37) 

Pgen given by eq. (!37|) is plotted in Figure [211 As expected the the sign of a determines whether 
the perturbation vertically shrinks or stretches the particle distribution. We also see that the system 
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\^N{N-l){N-2) 
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Figure 21: The unperturbed (solid) and perturbed (dotted) particle distributions, the latter given 
by eq. (1371) . for = 1000 and 9 = 1 with a) a = —j^, h) a = and c) a = 




Figure 22: The figures show the a) perturbed entanglement for a = 0.5 as a function of mo and 6 
for N = 100 and b) the difference - for = 100 and a = 0.5. 

is very sensitive to three-body recombination terms, as a coupling constant of order 0.001 causes 
significant changes to the particle distribution. 

Once again, setting P^^'^ equal to the last term in eq. (1371) . we can use eqs. (!24l) and (!25|) to 
study the perturbed entanglement. The choice of the sign of a to increase the entanglement again 
depends on the values of 9, Ai and A2. 

9 Discussion 

We have successfully studied the effects of a number of perturbations to the two-mode BEG model 
considered in [6]. We have found the corrections to the condensate wave functions, which in turn 
allowed the determination of the corrections to the particle distribution, time-evolution of the rel- 
ative number operator and the entanglement. In the non-degenerate case, we have shown that the 
model of [6], in which the coupling constants of the Hamiltonian are constrained, is robust to per- 
turbations in these constants. The system is most sensitive to perturbations in the elastic scattering 
length U and in the mode-exchange parameters A and fi. In each of the parameter perturbations 
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the entanglement of the coherent states |A^, N) and lA*", — A^) is most affected. Each parameter per- 
turbation was observed to increase the asymmetry in the particle distributions. Wc also observed 
that for specific values of mo and 6, the latter corresponding to certain parameter strengths, the 
entanglement perturbations are especially large. It was found that when the condensate is degen- 
erate (because of specific choices of Ai and A2) it is much more sensitive to perturbations, both in 
terms of particle distributions and entanglement. The effects on the entanglement are qualitatively 
different than in the non-degenerate case. In particular, the perturbations to the entanglement are 
mainly present only in the regions in which mo is close to —N which corresponds to a condensate 
with small negative scattering length. 

We have also extended the formalism to include the analysis of interactions involving particle 
loss. Prom these we can predict corrections to the particle distribution, entanglement and evolution 
of the relative population. This provides a new class of possible experiments that will allow the 
model here to be tested. Indeed, interactions involving particle loss have been limited experimen- 
tally, partially because they create instabihties in the system. With our results these interactions 
could be allowed to occur, and the results compared with the predictions contained above. The 
changes induced by both three-body recombination and background collisions are qualitatively dif- 
ferent than those induced by parameter perturbations. We find that the system is very sensitive 
to these external perturbations, which induce large changes in the particle distribution and en- 
tanglement from relatively small external coupling strengths, as compared to the induced changes 
from parameter perturbations. As with the degenerate parameter perturbation, we found that the 
perturbative effects on the entanglement become negligible as mo approaches A^, i.e. when the 
scattering length is small and positive. In general, we can conclude that higher collision rates make 
the condensate more stable to perturbations. 

Our results promise to be useful in the experimental realization of two-mode Bose-Einstein 
condensates which are stable to parameter perturbations and particle loss. We are planning to 
extend our analysis to include many-body interactions which are present at cooler stages of the 
condensate and study the role of such interactions in the stability of the condensate. 

Appendix: The Schwinger 5u(2) Boson Representation 

Let Jz and J± = J^, ± iJy be the usual generators of the Lie algebra su(2), satisfying 



[J„ J±] =±J±, [J+, J_] = 2J,. 




The eigenstates are labeled as \j, niang) where 
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Note that for fixed j, rUang may take any of the 2j + 1 values — j, —j + |, . . . , j — i, j. We also 
have J±\j,mang) = + 1) - rnangirriang ± l)|j,ma„g±l). The Schwinger representation of su(2) 

defines a Lie algebra homomorphism between the angular momentum representation (generated by 
J± and Jz) and a bosonic representation. Consider two bose operators a and h with [a, a^] = 1 = 
[6, 6^] and all other pairs having vanishing commutator. Defining the mapping 

J+ ^ a^6, J_ ^ah\ Jz ^ ^ {a^a - ¥b) (39) 

and extending linearly to the rest of su(2) then gives the desired homomorphism [21]. A short 
calculation then shows that is mapped to ^iV^ + where we have defined N = a^a + b^b. We 
label the basis states in the bosonic representation as | y, with the eigenvalue of N and m 
the eigenvalue of m = a^a — b^b; this is in complete analogy with the label \j,mang)- Note that 
f may take the 2 (f ) + 1 = + 1 values -f , -f + 1, . . . , f - 1, f . From the definition of 
the homomorphism it follows that a^b\f,f) = l^/N{N + 2) - m(m + 2) |f , f + l). If we then 
rescale the state | y, y) to \N,m) the above identity is rewritten as 

a^b\N,m) = ^v^A^(A^ + 2) - m(m + 2) |A^,m + 2) . (40) 

That I A^, m) is mapped to a multiple of |A^, m + 2) can be seen directly from the form of the operator 
a^b, which annihilates a particle in mode b while creating one in mode a. Similarly we have 

ab^ \N,m) = -^N{N + 2) — m{m — 2) |A^, m — 2) . (41) 

Now rescaled, m may take on the A^ + 1 values — A^, — A^ + 2, . . . , N — 2, N. 

Appendix: Counting Terms in the n-Model Hamiltonian 

It is of interest to quantify the generality of the n-model Hamiltonian under study. We begin 
this below by first counting the number of terms in the most general Hamiltonian that would be 
of interest to us. To do this we must define precisely the Hamiltonians that are of interest to 
our study of two- mode BECs. First, we limit ourselves to Hamiltonians that are polynomials in 
the bose operators. The interactions under consideration consist of any total number preserving 
operations. Operators corresponding to such interactions must then commute with the total number 
operator N. To restrict the class of interactions, we consider only those that have no intermediate 
interactions, such as the spontaneous creation and annihilation of a particle. This imposes the 
restriction that the Hamiltonian be normal ordered. As usual, we also require self-adjointness of 
the Hamiltonian. Finally, we may decompose the Hamiltonian into its homogeneous parts, i.e. 
terms of degree 1,2,3,.... So, it is sufficient to first consider only homogeneous Hamiltonians, 
and then construct more general Hamiltonians from these. With these assumptions we prove the 
following proposition. 
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Proposition 9.1. Let H be a homogeneous, self- adjoint, normal- ordered polynomial in the hose 
operators a and h and their adjoints. Furthermore, assume that each term in H commutes with the 
total number operator N = a^a + b^b. Put deg{H) = n > 0. Then the number of terms in H is at 
most (!^±Mn±4) . 

Proof. Let x(n) denote the number of terms in the Hamiltonian. Observe that the requirement that 
each term in H commute with ensures that n = 2k for some A; G N. Consider now a monomial of 
H with 2k — p of the operators being either a or while the remaining p being either 6 or 6^ It is 
easy to check that, assuming self-adjointness, normal ordering and vanishing commutator with A^, 
there are [|J + 1 such terms, where for g G M, \_q\ is the greatest integer less than or equal to q. We 
may group the monomials of H into three groups according to whether there are more than, less 
than, or the same number of a mode terms as b mode terms. Doing so, we find that the maximal 
number of terms in H is 

fc-i , 
x{n = 2k)=2j2{[l\+l) + [^\+l- 

Since [^j^J + [|J = k — l for all A; G Z we obtain the recursion relation x(2A;) = x(2(A; — 1)) + A; + 1. 
Repeated application of this relation yields x(2A;) = x(2) + XljisJ) "which we may rewrite as 

k+l 

x(2A;) = j since x(2) = 3, from which the proposition follows. □ 

Since we can decompose a general polynomial in terms of its monomials, the above proposition is 
sufficient to count the maximal number of terms in the most general Hamiltonian described above. 

Corollary 9.2. Let H be as above without the assumption of homogeneity with deg{H) = n > 0. 
Then the number of terms in H is at most -^n^ + + j^n + 1. 

Proof. Again, we can write n = 2k for some integer A; > 0. The number of terms in the Hamiltonian 
is the sum of the number of terms in each of its homogeneous parts, i.e. Yl'j=oX{'^j)- After some 
algebra the corollary follows. □ 

As noted in the proofs of the above proposition and corollary, the assumptions on the Hamilto- 
nian imply deg{H) = n is even. We see that the total number of terms in the most general n-model 
Hamiltonian grows like rt". We would like to compare this result to the number of terms in the 
n-model of [6l [7]. Unfortunately, we have not been able to successfully count the number of terms 
in the n-model. In order to do so, one would first need to conjugate 

n 

by the displacement operator U{^). The number of terms in this calculation grows quickly with 
n which makes the conjecture of a formula for the number of terms difficult. Moreover, general 
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arguments to count the number of terms, as in the proposition above, seem difficult to make in this 
case. Regardless, the results for n = 1,2 and 3, summarized in Table [H suggest that there are a 
significant number of terms missed by the n-model. Whether or not these terms are important in 
an experimental setting is another question. 

Table 1: Comparison of the number of terms in n-model and the most general model 



n 


Terms in ra-model 


Terms in general model 


Missed Terms 





1 


1 





1 


3 


4 


1 


2 


6 


10 


4 


3 


13 


20 


7 
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